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In  this  paper  a  new  multiobjective  modified  honey  bee  mating  optimization  (MHBMO)  algorithm  is 
presented  to  investigate  the  distribution  feeder  reconfiguration  (DFR)  problem  considering  renewable 
energy  sources  (RESs)  (photovoltaics,  fuel  cell  and  wind  energy)  connected  to  the  distribution  network. 
The  objective  functions  of  the  problem  to  be  minimized  are  the  electrical  active  power  losses,  the  voltage 
deviations,  the  total  electrical  energy  costs  and  the  total  emissions  of  RESs  and  substations.  During  the 
optimization  process,  the  proposed  algorithm  finds  a  set  of  non-dominated  (Pareto)  optimal  solutions 
which  are  stored  in  an  external  memory  called  repository.  Since  the  objective  functions  investigated  are 
not  the  same,  a  fuzzy  clustering  algorithm  is  utilized  to  handle  the  size  of  the  repository  in  the  specified 
limits.  Moreover,  a  fuzzy-based  decision  maker  is  adopted  to  select  the  ‘best’  compromised  solution 
among  the  non-dominated  optimal  solutions  of  multiobjective  optimization  problem.  In  order  to  see  the 
feasibility  and  effectiveness  of  the  proposed  algorithm,  two  standard  distribution  test  systems  are  used 
as  case  studies. 
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1.  Introduction 

The  application  of  the  RESs  such  as  wind,  fuel  cell  and  photo¬ 
voltaic  in  the  new  competitive  electric  power  markets  has  gained 
significant  attention  due  to  the  economic  and  environmental  con¬ 
cerns  of  fossils  and  nuclear  fuel-based  electricity  energy  as  well  as 
reduction  of  fossil  resources  [1  ].  Also,  the  existence  of  some  impor¬ 
tant  aspects  as  the  quality  of  the  RESs  such  as  compatibility  with 
other  modular  subsystem  packages,  fully  automation  possibility, 
low  emission  release,  high  efficiency  and  proper  power  quality  and 
reliability  have  made  them  even  more  popular  than  before  [2]. 

In  recent  years,  so  many  researchers  have  attended  to  investi¬ 
gate  the  use  of  some  kinds  of  renewable  energies  like  wind  energy, 
biogas  energy,  fuel  cells,  photovoltaic  cells,  combined  heat  and 
power  systems  (CHP),  etc.,  in  the  distribution  voltage  level  [3-6]. 
Nevertheless,  there  are  some  significant  considerations  to  get  use 
of  the  RESs  appropriately  and  efficiently.  Regions  like  offshore  and 
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high  altitude  areas  that  have  more  constant  and  stronger  winds 
are  suitable  to  be  used  for  the  construction  of  wind  farms.  The 
power  stored  in  the  airflows  can  be  employed  to  rotate  wind  tur¬ 
bines  and  so  generate  a  clean  and  consistent  electric  power.  Fuel 
cell  with  a  modular  structure  allows  for  simple  construction  and 
operation  with  possible  applications  for  distributed  and  portable 
power  generation  [7].  Also  as  a  result  of  their  fast  response,  fuel  cells 
have  a  good  quality  to  follow  and  supply  the  load  changes  while 
maintaining  the  high  efficiency  at  the  same  time  [3-6].  Another 
new  technology  in  the  field  of  renewable  energy  technologies  is 
photovoltaics  (PV).  PV  is  a  method  of  generating  electrical  power 
by  converting  solar  radiation  into  direct  current  electricity  using 
semiconductors  that  exhibit  the  photovoltaic  effect  [8].  Like  the 
other  kinds  of  renewable  energies,  PV  has  found  many  applica¬ 
tions  including  satellites,  electric  vehicles,  remote  dwelling,  boats, 
on  roofs,  and  by  the  use  of  DC-AC  converters  in  the  grids  which  are 
connected  to  the  power  system.  All  these  applications  and  many 
other  benefits  that  are  not  mentioned  here  make  it  critical  to  inves¬ 
tigate  the  effect  of  the  RESs  on  the  distribution  network  especially 
in  the  area  of  the  DFR  problem. 

Electric  distribution  networks  are  generally  designed  and  con¬ 
structed  as  the  radial  networks  so  as  to  have  suitable  and  proper 
protection  coordination.  Nevertheless,  the  necessity  of  having  a 
secure  network,  supplying  all  consumers,  minimizing  power  losses 
and  improving  power  quality,  it  is  required  to  change  the  struc¬ 
ture  and  the  topology  of  the  network  using  automatic  or  manual 


8882 


T.  Niknam  et  al.  /  Journal  of  Power  Sources  196(2011)  8881-8896 


Nomenclature 

X 

state  variables  vector 

n 

number  of  state  variables 

NFc 

number  of  FC  power  sources 

NpV 

number  of  PV  power  sources 

Nwind 

number  of  wind  power  sources 

Nb 

number  of  branches 

Ri 

resistance  of  zth  branch  (£2) 

h 

current  of  zth  branch  (A) 

Ppc  ,i 

active  power  production  of  the  zth  fuel  cell  power 
source  (kW) 

Pp  v,i 

active  power  production  of  the  zth  PV  power  source 
(kW) 

^Wind  ,i 

active  power  production  of  the  zth  wind  power 
source  (KW) 

^sub 

active  power  production  of  the  substation  (kW) 

m 

electrical  efficiency  of  the  zth  FC 

PLRj 

part  load  ratio  of  the  zth  FC 

Qc  ,i 

cost  of  electrical  energy  generated  by  of  the  zth  FC 
power  source  ($) 

Cpv.i 

cost  of  electrical  energy  generated  by  of  the  zth  PV 
power  source  ($) 

Qvind,/ 

cost  of  electrical  energy  generated  by  of  the  zth  Wind 
power  source  ($) 

Qub 

cost  of  power  generated  at  substation  bus  ($) 

Price 

cost  of  power  per  unit  generated  at  substation  bus 
($) 

Cr 

annual  rates  of  benefit 

LF 

loading  factor 

Efcj 

emission  of  the  zth  FC  power  source  (lb) 

Epv,i 

emission  of  the  zth  PV  power  source  (lb) 

^Wind.i 

emission  of  the  zth  wind  power  source  (lb) 

^Grid 

emission  of  large  scale  sources  (substation  bus  that 
connects  to  grid)  (lb) 

NOxFCiI- 

nitrogen  oxide  pollutants  of  the  zth  FC  power  source 
(lb  kWh-1) 

S02FC,i 

sulphur  oxide  pollutants  of  the  zth  FC  power  source 
(lb  kWh-1) 

NOxPV)2 

nitrogen  oxide  pollutants  of  the  zth  PV  power  source 
(lb  kWh-1) 

S02pv,i 

sulphur  oxide  pollutants  of  the  zth  PV  power  source 
(lb  kWh-1) 

NOxwind.i  nitrogen  oxide  pollutants  of  the  zth  wind  power 

source  (lb  kWh-1 ) 

so2Wind)i 

sulphur  oxide  pollutants  of  the  zth  wind  power 
source  (lb  kWh-1 ) 

NOxGrid 

nitrogen  oxide  pollutants  of  the  grid  (kg) 

S°2Grid 

sulphur  oxide  pollutants  of  the  grid  (kg) 

^min.FC,! 

minimum  active  power  of  the  ith  FC  power  source 
(kW) 

^rnax.FC  ,i 

maximum  active  power  of  the  ith  FC  power  source 
(kW) 

^Vnin.P  V,i 

minimum  active  power  of  the  zth  PV  power  source 
(kW) 

^max.PV.i 

maximum  active  power  of  the  zth  PV  power  source 
(kW) 

Pmin.wind,/  minimum  active  power  of  the  zth  wind  power 

source  (kW) 

Pmax.wind.i  maximum  active  power  of  the  zth  wind  power 

source  (kW) 

|pLme| 

absolute  power  flowing  over  distribution  lines  (kW) 

pLine 

1  ij,  max 

maximum  transmission  power  between  the  nodes  z 
and  j  (kW) 

pLine 

1  ij,  min 

minimum  transmission  power  between  the  nodes  z 

and  j  (kW) 

Vmax 

maximum  value  of  voltage  magnitudes  of  zth  bus  (V) 

^min 

minimum  value  of  voltage  magnitudes  of  zth  bus  (V) 

m 

zth  objective  function 

MX) 

equality  constraints  of  zth  objective  function 

a(x) 

inequality  constraints  of  zth  objective  function 

j^min 

lowest  limit  of  zth  objective  function 

jmax 

highest  limit  of  zth  objective  function 

Nf 

is  the  number  of  the  objective  functions  in  the  MOP 

lifiW 

membership  function  for  zth  objective  function 

D 

drone 

Xqueen 

best  particle  among  the  entire  population  or  the 
queen 

^broodj 

the  jth  brood 

Sp 

queen  spermatheca  matrix 

NsP 

size  of  the  queen  spermatheca 

m 

absolute  difference  between  the  fitness  of  the  drone 
and  the  fitness  of  the  queen 

a 

speed  reduction  factor 

y 

random  value  in  the  range  of  [0,1  ] 

Prob(D) 

probability  of  adding  the  sperm  of  drone  D  to  the 
queen  spermatheca 

SCO 

queen  speed 

FAX) 

values  of  the  augmented  /j(X) 

Ne  q 

number  of  equality  constraints  of  the  DFR  problem 

Nueq 

number  of  inequality  constraints  of  the  DFR  problem 

Li 

penalty  factor 

f  2 

penalty  factor 

^ipop 

number  of  the  bees 

Squeen 

queen  speed 

Smax 

maximum  speed  of  the  queen 

■^min 

minimum  speed  of  the  queen 

I<  i 

value  of  the  production  of  NOx  (lbk  Wh_1 ) 

K 2 

values  of  the  production  of  SOx  (lbK  Wh_1 ) 

rqueen 

•h 

the  value  of  the  zth  objective  function  for  the  queen 

^drone 

the  value  of  the  zth  objective  function  for  the  drone 

Wt 

the  weighting  of  the  zth  objective  function 

Mi 

the  mean  value  of  the  drones’  population  column¬ 
wise 

mt 

the  mean  value  of  the  zth  element  of  the  control 
vector  in  the  drones’  population  column-wise 

*k 

random  value  in  the  range  of  [0,1  ] 

Tf 

a  constant  factor  which  decides  the  value  of  mean 
to  be  changed.  Can  be  lor  2 

Xq,k 

the  kth  new  queen  generated  for  implementing 
modifying  the  breeding  process 

XD,m 

the  zrzth  new  drone  generated  for  implementing 
modifying  the  breeding  process 

round 

the  mathematic  function  which  rounds  each  value 
to  the  nearest  integer 

rand( ) 

the  function  for  the  generation  of  random  value 

Ykjn 

the  new  individual  generated  through  modification 
process 

Z 

the  new  individual  generated  through  modification 
process 

switches.  However,  the  radial  structure  of  the  networks  and  dis¬ 
crete  nature  of  the  switches  is  a  main  obstacle  to  get  use  of 
the  classical  optimization  methods  in  the  multiobjective  distribu¬ 
tion  feeder  reconfiguration  (MDFR)  problem.  Classical  optimization 
methods  have  suggested  transforming  the  multiobjective  opti- 
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mization  problem  to  a  single  objective  optimization  problem  [9,10]. 
However,  in  these  methods  extensive  computational  resources 
(memories)  are  needed  to  find  the  multiple  optima. 

In  the  area  of  MDFR  problem,  Das  [11]  has  presented  a  new  mul¬ 
tiobjective  approach  based  on  the  fuzzy  set  theory  to  solve  the  DFR 
problem.  Vanderson  Gomes  et  al.  have  suggested  a  new  heuris¬ 
tic  method  to  change  the  configuration  of  the  distribution  systems 
[12].  Kim  et  al.  [13]  have  suggested  a  new  method  based  on  neural 
network  in  order  to  identify  the  configuration  of  the  network  cor¬ 
responding  to  different  load  levels.  Taylor  and  Lubkeman  [7]  have 
presented  an  expert  system  to  reduce  the  search  space  by  the  use  of 
heuristic  rules.  Kashem  et  al.  [14]  have  introduced  a  new  algorithm 
based  on  the  distance  measurement  to  find  a  loop  in  the  network, 
and  then  a  switching  plan  is  used  to  enhance  the  load  balancing  in 
the  mentioned  loop.  Baran  and  Wu  [15]  have  used  a  mixed  integer 
programming  (MIP)  to  solve  the  load  balancing  problem  and  power 
loss  reduction  in  the  network.  In  order  to  supply  the  purpose  of  load 
balancing  and  service  restoration  in  two  feeder  networks,  an  algo¬ 
rithm  based  on  the  network  reconfiguration  has  been  presented  by 
Zhou  et  al.  [16].  Also,  in  [17,18]  Niknam  has  proposed  a  method 
based  on  the  evolutionary  algorithms  to  convert  multiobjective 
DFR  problem  to  an  equivalent  augmented  single  objective  prob¬ 
lem.  Here  again  too  extensive  computational  resources  and  many 
runs  are  needed  to  find  the  solutions;  therefore  the  efficiency  of  the 
algorithm  will  be  reduced.  Also,  the  new  hybrid  algorithm  used  by 
Niknam  in  [19]  to  solve  the  multiobjective  DFR  problem  results  in  a 
single  solution  which  does  not  observe  the  Pareto  optimality  con¬ 
cept.  In  fact  it  has  the  main  disadvantage  of  ignoring  many  other 
good  candidate  solutions  that  can  be  supposed  as  optima. 

In  this  paper  a  new  modified  honey  bee  mating  evolutionary 
algorithm  is  utilized  to  solve  MDFR  problem  while  the  effect  of 
RESs  is  considered.  The  traditional  HBMO  suffered  from  two  main 
shortcomings:  (I)  dependency  on  the  algorithm  parameters,  and  (II) 
the  possibility  of  being  trapped  in  local  optima.  Therefore,  in  this 
paper  a  new  modification  process  is  proposed  to  enhance  the  per¬ 
formance  of  the  algorithm.  Also,  during  the  optimization  process, 
the  set  of  obtained  non-dominated  solutions,  called  Pareto-optimal 
solutions,  are  stored  in  the  repository.  In  order  to  control  the  size 
of  the  repository,  a  fuzzy  clustering  technique  is  utilized. 

2.  Multiobjective  DFR  considering  RESs 

2.1.  Objective  functions 

-  Minimization  of  the  power  losses :  Total  power  losses  can  be  mini¬ 
mized  by  the  following  equation: 

Ar 

MX)  =  P,„SS(X)  =  ^ R<  x\'i\2’X  =  [Tie.  Sw, P„] 

Tie  =  [Tie, ,  Tie2 ,  Tie3 , . . . ,  TieNtie  ] ;  Sw  =  [Sw, ,  Sw2 ,  Sw3 , . . . ,  SwNsw  ]  ( ? ) 

Ac  —  [Ac,i  ,  PfC,2  ,  •  •  • ,  Pfc,n¥C  ] :  Ppv  —  [  Av,l ,  Ppv,2 ,  .  •  •  ,  PpV,NpV  ] 

Pg  =  [Ac,  Ppv,  A/Vind] 

Pwind  =  [A/Vind,l  >  A/Vind,2>  •  •  •  >  A/Vind,Nwind 

where  Tie/  and  Sw*  are  the  states  of  the  ith  tie  switch  and  sec- 
tionalizing  switch  which  0  and  1  are  the  values  corresponding  to 
open  and  close  states,  respectively. 

-  Minimization  of  the  voltage  deviation  of  the  buses :  This  objective 
function  can  be  defined  as  follows: 

/2(X)  =  dev(X)  =  max[|l  -Vmin|  and  |l-Vmax|]  (2) 

Minimization  of  the  total  cost  of  generation:  The  total  cost  is  the 
summation  of  the  cost  related  to  the  power  produced  by  the  grid 
and  the  cost  related  to  the  power  produced  by  the  RESs.  The  grid 
cost  can  be  evaluated  as  follows: 


Table  1 

Specifications  of  RESs  (Test  system  1). 


Capacity  (kW) 

Type 

Location 

RES  1 

464.375 

FC 

6 

RES  2 

464.375 

FC 

13 

RES  3 

464.375 

PV 

19 

RES  4 

464.375 

Wind  energy 

22 

The  cost  of  FC  power  sources  can  be  evaluated  as  follows  [20]: 
CFC  i  =  0.04 $kWh-1  x  — 

m 


1  maXj 

if  PLR,  <  0.05  =>  iji  =  0.2716 

if  PLR,  >  0.05  =x  t]t  =  0.9033  PLR,5  -  2.9996 PLRf  +  3.6503  PLR? 
-2.0704  PLR?  +0.3747 

The  cost  of  PV  and  wind  units  can  be  evaluated  similarly  by  Eq.  (5). 
The  cost  of  generation  of  each  kWh  is  a  function  of  three  parameters 
[21]:  (I)  investment  cost;  (II)  operation  and  maintenance  cost;  (III) 
fuel  cost. 

Cpv,i  =  a  +  b  x  Ppv;f 
Q/Vind,i  =  a  -f  b  x  Pwind,  i 

Capital  cost($kW-1 )  *  Capacity(kW)  *  Gr  (5) 

a=  Life  time(Year)  *  365  *  24  *  LF 

b  =  Fuel  cost($kWh-1 )  +  0&MCost($kWh_1 ) 

Therefore  the  total  cost  is  as  follows: 

NFC  Npv  Nwind 

/3(X)  =  Cost  =  5>,  +  ]Tcpv,;  +  Cwind,i  +  csub  (6) 

i=l  1=1  i=l 


Minimizing  the  total  emission  produced :  The  total  emission  of  the 
grid  and  the  RESs  is  as  follows: 

nfc  Npv  Nwind 

/4PO  =  Emission  =  ^~^EfC,/  +  V,i  +  ^Wind ,i  +  ^Grid 


i=l  i=l  i=l 

Fr  i  F c  i  lb  MWh-1 

Efc,,  =  NOXpc,,  +  S02fc,,  =  (X[CJ  +  Kf )  x  PFC, 

PV  /  pv  1  lb  MWh  l 

E pV,/  =  NOxpv.f  +  S02PV  j  =  {K™>1  +  K™’1)  x  PPV,i 

c  1  cm  /r^Wind,i  rewind, id*3  MWb 

^Wind,i  —  In  Unwind,  i  +  oUZwind,i  —  +  ^2  ' 

x^Wind,i 

Ec rid  =  NOxCrid  +  S02Grid  =  (X“d  +  K2Grid)'b  MWh_’  x  Psub 


(7) 


The  relevant  values  of  these  parameters  are  brought  in  Table  2. 


2.2.  Limits  and  constraints 


(i)  Limits  associated  with  distribution  lines: 

nLine  ^  pLine  <  pLine 
1  y,min  ^  1  ij  ^  1  y,max 

(ii)  Distribution  power  flow  equations: 

Nbus 

Pi  =  JfViVjYij  COS (%-5,+Sj) 

1=1 

Nbus 

Q,  =  J2* (i) (ii) * * viyjYij  sin (ffy-St  +  Sj) 

1=1 


(8) 


(9) 


This  equation  is  a  load  flow  equation  that  plays  the  role  of  an  equal¬ 
ity  constraint. 


^sub  —  price  x  P$ub 


(iii)  Keeping  the  radial  structure:  since  the  distribution  networks 
(3)  are  supposed  to  be  radial,  this  quality  of  the  network  should  be 
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Table  2 

Emission  factors  related  to  N0X,  C02  and  S02. 


Emission  factors  (lb  MWh-1 ) 


Emission  type 

Grid 

Gas  turbine 

Micro  turbine 

Wind 

PV 

FC 

Internal  combustion  (IC) 

NOx 

5.06 

0.03 

0.44 

0 

0 

1.15 

4.7 

co2 

2031 

1078 

1596 

0 

0 

1108 

1432 

so2 

7.9 

0.006 

0.008 

0 

0 

0.008 

0.454 

preserved  during  the  reconfiguration.  Each  loop  in  the  network 
is  composed  of  a  tie  switch  and  a  sectionalizing  switch  that 
each  time  that  one  is  open;  the  other  one  is  closed,  so  that  a 
radial  network  would  be  achieved. 

(iv)  Feeder  current  limitation  [22]: 

\Ifi\<I™x>  i  =  h2,...,Nf  (10) 

(v)  RESs  constraints  on  active  power  production: 

Pmin,FC,i  —  PFC,i  —  Pmax,FC,i  *  Pmin,PV,z  —  PPV,i  —  Pmax,PV,z* 

Pmin,Wind,z  —  Pwind,i  —  Pmax,Wind,i  (H) 

(vi)  Bus  voltage  constraints: 

Vmin  <  V*  <  Vmax  (12) 

3.  Modeling  of  the  RESs 

3.1.  Fuel  cell 

FC  has  been  among  of  the  most  important  developments  in  the 
field  of  power  generation  in  recent  years.  Because  of  their  clean¬ 
ness,  good  efficiency  and  high  reliability,  they  have  become  one  of 
the  most  attractive  power  supply  sources  in  the  field  of  distributed 
generation  and  electrical  vehicles  [23,24].  There  are  many  different 
kinds  of  FCs  according  to  their  different  characteristics.  The  solid 
oxide  fuel  cells  (SOFCs)  and  molten  carbonate  fuel  cells  (MCFCs)  are 
supposed  to  be  in  the  high-temperature  category  of  the  FCs.  Their 
higher  efficiency  and  lower  operation  cost  is  considerable  in  com¬ 
parison  to  the  conventional  power  plants  in  the  megawatt  range  of 
production  [24].  The  typical  efficiency  of  conventional  power  plants 
ranges  from  38%  to  40%,  whereas  the  efficiency  of  an  SOFC  is  in  the 
range  of  55-60%  [25].  Another  kind  of  FC  is  the  proton  exchange 
membrane  FC  (PEMFC)  which  has  the  ability  of  operation  at  the  air 
temperature  allowing  rapid  startup  [23].  The  PEMFC  efficiency  is 
around  40-60%  with  the  power  and  voltage  ability  to  meet  special 
demands  as  a  result  of  its  modularity  structure. 

3.2.  Photovoltaic 

Solar  photovoltaics  (PVs)  are  arrays  of  cells  containing  a  material 
that  converts  solar  radiation  into  direct  current  (DC)  electricity  [8]. 
Solar  photovoltaic  (SPY)  electric  power  generation  is  a  promising 


clean  technology  with  vast  potential  [26].  Orbiting  satellites  were 
the  first  practical  application  of  PVs.  But  today  so  much  attention 
is  paid  to  use  them  as  a  source  of  power  in  the  power  generation 
grids.  As  the  PV  cell  temperature  exceeds  a  threshold  of  45  °C,  the  PV 
performance  decreases  [8].  In  this  situation  the  PVs  performance  as 
renewable  energy  power  source  should  be  investigated  in  different 
condition  carefully. 

3.3.  Wind 

The  performance  of  the  wind  turbines  is  based  on  the  conversion 
of  the  kinetic  energy  of  wind  into  electricity.  The  instantaneous 
power  produced  by  a  turbine  is  proportional  to  the  third  power 
of  the  instantaneous  wind  speed  [8].  Therefore  as  the  wind  speed 
increases,  the  output  power  increases  dramatically.  As  mentioned 
before,  regions  like  offshore,  open  flat  areas  and  high  altitude  areas 
that  have  more  constant  and  stronger  winds  are  suitable  to  be  used 
for  the  construction  of  wind  farms. 

3.4.  Modeling  of  the  renewable  energy  power  plants 

RESs  can  be  modeled  by  two  types:  (I)  constant  voltage  magni¬ 
tude,  constant  active  power  loads;  (II)  constant  active  and  reactive 
power  loads.  When  the  renewable  energy  power  source  is  supposed 
to  be  constant  active  power  and  voltage  magnitude  load  model,  it 
should  be  able  to  keep  its  voltage  magnitude  constant  by  the  injec¬ 
tion  of  reactive  power.  When  the  renewable  energy  power  source 
is  considered  as  a  constant  active  and  reactive  power  load  model,  a 
specific  active  and  reactive  power  is  produced  and  injected  to  the 
network.  However,  in  both  cases  it  should  be  considered  that  active 
and  reactive  power  produced  should  not  exceed  the  generation 
capacity.  As  the  result  of  the  structure  of  the  distribution  networks, 
the  load  distribution  in  these  networks  is  unbalanced.  So  the  RESs 
operation  and  control  is  done  in  two  forms:  (I)  simultaneous  three 
phase;  (II)  independent  three-phase  control  or  single  phase  con¬ 
trol.  Therefore  according  to  their  model  and  type  of  control,  four 
models  can  be  defined  for  these  generators  (Fig.  1)  [7]: 

(I)  Constant  active  and  reactive  load  model  with  simultaneous 
three-phase  control. 

(II)  Constant  reactive  power  and  voltage  magnitude  load  model 
with  simultaneous  three-phase  control. 


(b) 


6 


a 


6C 


Fig.  1.  Models  of  RESs:  (a)  constant  active  and  reactive  load  model  with  simultaneous  three-phase  control,  (b)  constant  active  and  reactive  load  model  with  independent 
three-phase  control,  (c)  constant  active  power  and  voltage  magnitude  load  model  with  simultaneous  three-phase  control  and  (d)  constant  active  power  and  voltage  magnitude 
load  model  with  independent  three-phase  control. 
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Fig.  2.  Block  diagram  of  MHBMO  algorithm. 
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Table  3 

Comparison  of  active  power  losses  objective  functions  evaluated  by  different  methods  neglecting  RESs  (Test  system  1 ). 


Method 

Power  loss  [kW] 

Minimum  voltage 

Open  switches 

Goswami  and  Basu  [29] 

139.53 

0.93781964 

57,59,514,532,537 

Vanderson  Gomes  et  al.  [30] 

139.53 

0.93781964 

s7,s9,sl4,s32,s37 

PSO-SFLA  [31] 

139.53 

0.93781964 

s7,s9,sl4,s32,s37 

MSFLA  [32] 

139.53 

0.93781964 

s7,s9,sl4,s32,s37 

Shirmohammadi  and  Hong  [33] 

140.26 

0.93781964 

s7,sl0,sl4,s32,s37 

The  proposed  algorithm 

139.53 

0.93781964 

s7,s9,s!4,s32,s37 

Table  4 

Comparison  of  voltage  deviation  objective  functions  evaluated  by  different  methods  neglecting  RESs  (Test  system  1). 

Method 

Voltage  deviation  [p.u] 

Minimum  voltage 

Open  switches 

GA 

0.06218097 

0.93781902 

s7,sl0,sl4,s32,s37 

PSO 

0.06120031 

0.93879681 

s6,s9,sl4,s32,s37 

HBMO 

0.06266643 

0.93733335 

s9,sl2,sl9,s35,s37 

The  proposed  algorithm 

0.06120031 

0.93879681 

s7,s9,s!4,s32,s37 

(III)  Constant  active  and  reactive  load  model  with  independent 
three-phase  control. 

(IV)  Constant  reactive  power  and  voltage  magnitude  load  model 
with  independent  three-phase  control. 


such  that  X  dominates  X*  e  £2.  £2  is  the  set  of  all  the  vectors  (X)  that 
observe  the  constraints  and  limitations. 

In  definition  the  solution  X\  dominates  X2  if  the  following  two 
conditions  are  satisfied: 


4.  Methodology 

Multiobjective  optimization  problem  (MOP)  is  the  process  of 
optimizing  some  different  conflicting  objective  functions  simulta¬ 
neously  when  all  the  constraints  and  the  limitations  are  observed 
and  the  best  optimized  solution  for  all  the  objective  functions  are 
achieved.  MOP  can  be  defined  as  [22]: 

minF  =  Lfi(X),/2(X),...,/„(X)]r 

s.t.  (13) 

fi(X)<0  1  =  1,2,  ...,Nueq  1  ' 

h,(X)  =  0  1  =  1,2,  ...,Neq 

In  this  essay  X  is  the  variable  vector  of  making  decision.  Also  n  is 
the  number  of  objective  functions. 

X  =  [Tie,Sw,Pg]  (14) 

In  fact  the  main  difference  between  the  single  and  the  multiobjec¬ 
tive  optimization  problem  is  the  ability  of  the  MOP  in  observing 
the  best  optimal  solution  of  different  objective  functions  simulta¬ 
neously.  This  feature  is  due  to  capability  of  MOP  in  selection  of  the 
solutions  as  a  set  of  Pareto  points.  In  MOP  problem  X*  is  called  a 
Pareto  optimal  solution  if  it  is  impossible  to  find  a  solution  X  in  £2 


(1  )Vje  {l,2,...,n},J5(X,)<J5(X2) 
(2)3  k  e  {1,2  ,...,n},/k(X,)</k(X2) 


5.  Fuzzy-based  clustering 


When  the  optimization  is  started,  in  each  iteration,  the  best 
values  of  the  objective  functions  are  evaluated.  However,  the  big 
difference  which  exists  between  the  optimal  value  of  some  of  the 
objective  functions  and  the  optimal  value  of  the  other  ones  is  a 
barrier  for  making  the  improvement  rate  similar  to  each  other.  The 
main  point  of  using  fuzzy  set  theory  in  this  investigation  is  to  bring 
all  the  optimal  values  of  the  objective  functions  in  the  same  base, 
so  that  to  create  a  good  criterion  for  comparison.  The  membership 
function  designated  to  each  objective  function  is  as  follows: 


1 

0 


fmax 

H _ 

fmax 

H 


-/i(X) 

_  fmin 


for  MX)<f™* 
for  /i(X)>^max 

frn  </«(x)  </;max 


(16) 


As  we  see  the  membership  functions  used  here  are  continuous 
functions  with  lower  and  upper  limits  which  decrease  monoton- 
ically. 


Table  5 

Comparison  of  objective  functions  evaluated  by  different  methods  considering  RESs  for  25  trails  on  Test  system  1. 


Objective  function 

Method 

Average 

Standard  deviation 

Worst  solution 

Best  solution 

Open  switches  of  best  solution 

CPU  time  [s] 

Power  losses  [kW] 

GA 

95.435466 

3.774 

104.12418 

92.273284 

s7,sl0,s31,s34,s37 

13.924 

PSO 

89.975333 

1.916 

95.265620 

89.089835 

s6,s8,s32,s34,s37 

12.785 

HBMO 

92.548024 

2.115 

98.333334 

91.103298 

s6,s8,s34,s36,s37 

13.234 

The  proposed  algorithm 

85.583101 

00.00 

85.583101 

85.583101 

s7,sll,s31,s34,s37 

10.031 

Voltage  deviation 

GA 

0.05487009 

0.002391 

0.06104787 

0.05291578 

s7,s9,sl3,s32,s37 

11.593 

[p.u] 

PSO 

0.05083480 

0.001547 

0.05444417 

0.04920617 

sll,sl3,s35,s36,s37 

10.319 

HBMO 

0.05120344 

0.0021054 

0.05851443 

0.05114542 

s6,sl3,s21,s32,s37 

11.043 

The  proposed  algorithm 

0.04898826 

0.0000000 

0.04898826 

0.04898826 

si  1  ,s31  ,s33,s34,s37 

9.303 

Cost  [$] 

GA 

155.16237 

0.426 

155.85994 

154.98564 

s6,s9,s34,s36,s37 

12.776 

PSO 

154.43350 

0.235 

155.07246 

154.33555 

s7,sll,sl4,s36,s37 

11.350 

HBMO 

154.68270 

0.275 

155.35477 

154.39087 

s7,sll,s34,s36,s37 

11.945 

The  proposed  algorithm 

154.18182 

0.000 

154.18182 

154.18182 

s7,s9,sl4,s32,s37 

10.500 

Emission  [lb] 

GA 

27,587.566 

362.34 

28,046.59 

27,298.962 

s21,s28,s33,s32,s37 

12.245 

PSO 

27,376.903 

347.63 

26,878.888 

25,747.949 

sl2,s20,s35,s36,s37 

10.694 

HBMO 

27,440.376 

351.38 

27,523.738 

26,357.004 

s9,s33,s34,s36,s37 

11.404 

The  proposed  algorithm 

25,280.766 

000.00 

25,280.766 

25,280.766 

s7,s8,s32,s35,s37 

9.484 
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Fig.  3.  Single  line  diagram  of  32  bus  test  system. 


For  each  of  the  solutions  in  the  repository,  the  normalized  mem¬ 
bership  function  can  be  evaluated  as  follows: 


N„U)  = 


Eli x 

5Zi=l  x  ) 


(17) 


where  n  is  the  number  of  the  objective  functions  and  m  is  the  num¬ 
ber  of  the  individuals  in  the  repository.  This  equation  makes  a  type 
of  decision  making  criteria  which  is  adaptive  and  has  the  ability 
of  applying  the  decision  maker’s  options.  For  the  controlling  the 
size  of  the  repository,  all  the  weighting  factors  (w/)  are  set  unit. 
Therefore  for  all  the  non-dominated  solutions  is  evaluated  and 
the  best  solutions  are  kept  in  the  repository.  In  this  paper  we  get 
advantage  of  fuzzy  clustering  in  another  application  too.  In  the 
simulation  results,  it  is  explained  that  by  changing  the  weighting 
factors  to  a  different  value  than  1,  the  preferences  of  the  operator 
in  determining  the  importance  of  the  objective  functions  can  be 
achieved. 


generating  the  new  broods,  they  should  be  protected  and  improved 
by  the  workers  (for  a  complete  description  see  [27,28]).  As  a  matter 
of  fact,  the  new  generation  is  consisted  of  all  the  new  solutions  that 
should  be  used  in  the  optimization  process.  However,  if  one  of  these 
broods  has  a  better  situation  than  the  queen’s,  then  it  replaces  the 
queen  and  another  new  generation  is  created  by  the  use  of  the  new 
queen.  This  process  is  continued  until  the  time  that  the  best  queen 
is  achieved. 

6.2.  Modified  HBMO  (MHBMO)  algorithm 

As  we  mentioned  earlier,  original  HBMO  suffers  from  the  proba¬ 
bility  of  being  trapped  in  local  optima.  This  shortcoming  roots  from 
the  mating  process.  As  we  said  in  the  last  section,  the  process  of 
mating  between  a  drone  and  a  queen  in  the  original  HBMO  is  done 
probabilistically  by  using  Eq.  (18).  However,  in  the  MHBMO  algo¬ 
rithm,  each  queen  value  is  composed  of  the  all  value  of  the  objective 


6.  Honey  bee  mating  optimization  (HBMO)  modeling 

6.1.  Original  HBMO 

In  the  natural  world,  these  types  of  flying  insects  (honey  bees) 
live  with  each  other  as  a  colony.  That  is  their  life  has  a  direct  rela¬ 
tion  with  their  social  colony.  The  communication  of  these  insects 
structurally  is  composed  of  three  main  groups:  the  queen  (female), 
the  drones  (male)  and  the  workers.  The  process  of  mating  between 
a  drone  and  a  queen  is  done  probabilistically  by  using  an  annealing 
function  as  follows  [27]: 


If  the  mating  process  is  done  successfully,  the  drone  sperm  is  stored 
in  the  queen  spermatheca.  After  each  mating,  the  queen  speed 
decreases  as  follows: 

S(t  +  l)  =  axS(t)  (19) 

This  mating  process  continues  until  the  time  that  the  queen’s  sper¬ 
matheca  becomes  full  or  her  speed  deceases  to  a  specific  value.  After 


Table  6 

Some  of  the  non-dominated-solution  found  for  MDFR  problem  (Test  system  1 ). 


Power  losses  (kW) 

Voltage  deviation  (p.u) 

Cost  ($) 

Emission  (lb) 

1 

085.583 

0.053385 

158.060 

26,215.79 

2 

092.296 

0.054754 

157.705 

29,083.56 

3 

089.471 

0.058755 

157.026 

31,569.95 

4 

119.517 

0.060206 

155.648 

36,982.22 

5 

108.316 

0.055584 

156.445 

34,074.18 

6 

094.805 

0.048988 

158.429 

25,335.31 

7 

141.097 

0.066092 

154.437 

48,966.26 

8 

106.231 

0.055065 

156.207 

35,622.73 

9 

097.291 

0.049462 

157.391 

31,977.56 

10 

094.752 

0.056561 

157.840 

27,501.51 

11 

121.030 

0.062998 

159.305 

26,457.16 

12 

113.175 

0.056907 

155.816 

35,854.27 

13 

098.836 

0.057267 

156.827 

33,453.47 

14 

139.533 

0.062180 

154.181 

49,954.62 

15 

100.242 

0.055331 

158.113 

27,610.90 

16 

103.064 

0.053654 

156.429 

34,186.32 

17 

090.596 

0.054528 

158.260 

25,280.76 

18 

108.522 

0.060597 

156.552 

32,554.66 

19 

124.000 

0.057638 

155.324 

38,607.39 

20 

112.660 

0.056966 

155.640 

38,336.96 
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functions  which  are  considered  in  the  problem.  Therefore  the  in  the 
proposed  method,  Eq.  (18)  is  developed  as  follows: 


prob(D)  =  exp  (  - 

A/=  Vtf? 


m 


queen  _  f drone  _|_  ^-queen  _ydrone^2  _|_  ^-queen  _ydrone^2  _|_ 


(20) 


y drone  y 


In  the  original  HBMO  each  brood  is  produced  by  the  following  equa¬ 
tion: 


Xqueen  —  [*q,  1 ,  2,  •  •  •  ?  *q,N] 

SPi  =  [su,si)2,  ...sUN],i  =  1,2,  ...,NSp  (21) 

^broodj  —  Xqueen  +  Y  x  (Xqueen  —  $Pi) 

As  we  can  see,  in  the  original  HBMO,  one  queen  will  mate  with  the 
drones’  population  to  generate  the  new  broods.  But  in  the  proposed 
MHBMO  algorithm,  in  order  to  enhance  the  diversity  of  the  new 
generation  of  the  honey  bees  (broods),  after  generation  of  the  queen 
spermatheca  in  a  similar  method  to  the  original  HBMO,  the  breed¬ 
ing  process  is  corrected  as  the  following.  Note  it  that  this  process 
should  be  repeated  for  all  the  drones  (Xj). 

Firstly,  three  new  individual  should  be  generated  as  the  new 
queens  for  implementing  breeding  process.  The  first  individual 
(Xqi )  is  selected  from  the  repository  randomly.  For  the  generation 
of  the  second  individual,  the  mean  value  of  the  repository  is  calcu¬ 
lated  column-wise,  which  gives  the  mean  value  of  each  particular 
individual  as: 


Mf  =  [mum2,m3,  ...,mn]  (22) 

Note  it  that  Mj  which  is  generated  here  is  selected  as  the  new  sec¬ 
ond  queen  (Xq2).  The  third  modified  queen  (Xq3)  is  generated  by 
the  use  of  fuzzy  membership  function.  That  is  the  individual  in  the 
repository  which  the  summation  of  its  membership  function  val¬ 
ues  is  the  most  (so  the  best  solution)  is  the  selected  as  the  third 
modified  queen. 

After  selection  of  the  new  three  queens,  two  individual  should  be 
generated  as  the  new  drones  for  implementing  breeding  process. 
The  first  new  drone  is  the  mean  value  of  the  queen  spermatheca 
which  is  calculated  column-wise.  The  second  new  drone  is  selected 
from  the  population  of  the  drones  randomly.  The  idea  which  is  used 
here  for  the  modification  process  is  to  move  Xj  toward  the  corre¬ 
sponding  state  vector  of  the  queen.  In  fact  since  the  queen  position 
is  the  best,  so  by  moving  the  position  of  each  individual  toward 
the  queen’s  position,  the  position  of  the  specified  individual  will  be 
improved.  Now,  by  the  use  Xqi ,  Xq2,  Xq3,  Xm  and  XD2,  six  AX,-(  AXq  i-, 
A XD,i)  are  produced  as  the  follows: 

AXqk  Dm  —  ^k(^q,k  ~  TpXp)  m),  k  =  1, 2,  3,  171  =  1,2  (23) 

Here  TF  is  a  constant  factor  which  decides  the  value  of  mean  to  be 
changed.  Tp  can  takes  two  values  of  1  and  2  which  is  determined 
heuristically  as:  Tp  =  round[l  +rand(0,l)].  For  improving  the  diver¬ 
sity  of  the  search  space  a  drone  Xj  is  selected  from  the  population 
of  the  drones  such  that  Xj  =f=  X,-.  Then  by  the  use  of  the  following 
equation  A  Xj  is  determined  as: 

if  f(Xj)  >/(Xj) 

AXj=rand(.)x(XJ-XI) 

else  (24) 

A  Xj  =  rand(.)  x  (Xj  -Xj) 
end 


where  rand(.)  is  a  random  function  generator.  The  new  six  modified 
individuals  are  generated  as  follows: 

YKm=Xfd  +  AXqKDm  ■  k=  1,2,3,  m  =  1, 2 
Z  =  X°ld  +  AXj  1  j 

Now  by  the  use  of  Eq.  (15),  the  non-dominated  solutions  among  Y\t\ , 
Yi)2»  Y2ii  ,  Y2  2,  Y3ii  ,  Y3ii  and  Z  are  evaluated  and  stored  in  the  repos¬ 
itory.  If  there  is  only  one  non-dominated  solution  among  them, 


then  this  individual  is  selected  as  the  modified  brood.  Neverthe¬ 
less  if  there  is  more  than  one  non-dominated  solution  among  the 
individuals,  then  the  individual  which  the  summation  of  its  mem¬ 
bership  function  is  the  most  (so  the  most  satisfying)  is  selected  as 
the  new  modified  brood.  For  every  individual  in  the  population  of 
drones,  this  process  should  be  repeated.  After  that  the  breeding  pro¬ 
cess  for  all  the  drones  is  completed  then  the  repository  is  updated 
by  Eq.  (15). 

Another  significant  difference  between  the  MHBMO  algorithm 
and  the  original  HBMO  algorithm  is  in  the  generation  of  the  new 
drones’  population.  In  the  original  HBMO,  after  choosing  the  new 
queen  from  the  broods’  population  then  the  old  drone  generation 
is  discarded.  In  the  MHBMO  algorithm,  each  individual  of  the  new 
drones’  population  is  generated  during  the  breeding  process.  As 
mentioned  before,  for  each  drone,  seven  new  modified  broods  are 
generated.  After  each  breeding  process,  the  non-dominated  solu¬ 
tion  which  the  summation  of  its  membership  function  values  is 
the  most  (and  therefore  the  best  brood)  is  compared  with  that  of 
the  corresponding  drone  (Xj).  If  the  summation  of  the  best  brood 
membership  function  is  better  than  that  of  Xit  then  replace  it,  else 
Xj  will  be  kept  in  its  position.  So  after  a  complete  breeding  process, 
the  new  drones’  population  is  updated.  The  MHBMO  algorithm  is 
depicted  in  Fig.  2. 

6.3.  Solution  procedure  ofMDFR  considering  RESs 
Defining  the  input  data. 

Converting  the  constrained  MOP  to  an  unconstrained  one: 
here  the  constrained  MOP  should  be  changed  to  an  uncon¬ 
strained  one  by  constructing  an  augmented  objective 
function  as  Eq.  (26).  Since  all  constrains  must  be  met,  so 
penalty  factors  LI  and  L2  are  used  to  prevent  violating  the 
constraints. 

'Fl(X)' 

F2(X) 

F3(X) 

.  M  J  4x1 

Neq 

MX)  +  l^W))2  +  l2 
1=1 
Neq 

/2(X)  +  L1^a!(X))2+L2 
1=1 
Neq 

/3(X)  +  L1^ai(X))2+L2 
1=1 
Neq 

/4(x)  +  L1^yi(x))2  +  t2 

1=1 

(26) 

L\  and  L2  as  the  penalty  factors  are  supposed  to  be  105  in  the  paper. 

Step  3:  Initial  population  generation.  The  initial  population  is  as 
follows: 


^(Max[0,  -gj(X)])2 
1=1 
Nueq 

^(Max[0,  -g,(X)])2 
1=1 
Nueq 

^(Max[0,  -g,(X)])2 
1=1 
Nueq 

^(Max[0,  -gj(X)])2 


Step  1 : 
Step  2: 


F(X)  = 


initial-population  = 


Xi 

X2 


XN 


Njpop  x  (Nqe  +Nsw+Ng ) 


(27) 


Xj  =  [Xj]  =  [Tie;.  Swj,  PgA\ 
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Fig.  4.  3D  plot  of  the  Pareto-optimal  solutions  found  for  case  study  1  considering  total  active  power  losses,  voltage  deviation  and  cost  objective  functions. 


Tie  =  [Tiej, Tie2, Tie3, . . . ,  Tie/vTie], 

Sw  =  [Swi,  Sw2,  Sw3, . . . ,  SwjvSw] 


^wind,NFC  L  i  =  1 »  2,  3, . . . ,  N;  Ng  =  NPV  +  NFC  +  Nwind 


Pg  =  [^FC,  ^PV,  ^WindL  PFC  =  [^FC,1 ,  ^FC,2,  •  •  •  ,  PfC,Nfc  ] 

PpV  =  [Ppv,\ ,  ^PV,2,  •  •  •  ,  PpV,NPVl  Pwind  =  l^Wind,!  >  ^Wind,2> 


Step  4:  Objective  function  evaluation.  After  generating  the  initial 
population,  load  flow  is  carried  on,  and  all  the  membership 
functions  related  to  the  objective  functions  are  evaluated 
by  Eq.  (16). 


Fig.  5.  3D  plot  of  the  Pareto-optimal  solutions  found  for  case  study  1  considering  total  active  power  losses,  voltage  deviation  and  emission  objective  functions. 
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Fig.  6.  3D  plot  of  the  Pareto-optimal  solutions  found  for  case  study  1  considering  total  active  power  losses,  cost  and  emission  objective  functions. 


Here  /Xploss  j(X),  Mdvoitage,iP0>  Mprice,iPO  /^Emission,! P0 

membership  functions  of  the  total  power  losses,  the  voltage  devi¬ 
ation,  the  total  cost  and  the  total  emission  for  the  ith  individual  in 
the  population. 


Step  5:  Repository  complement.  Here  by  the  use  of  membership 
functions  evaluated  in  the  last  step  and  by  the  use  of  Eq. 
(15),  all  the  non-dominated  solutions  are  found  and  stored 
in  the  repository. 

Step  6:  Queen  selection.  A  queen  is  selected  from  the  non- 
dominated  solutions  stored  in  the  repository  randomly. 


Step  7:  Generation  of  the  queen  spermatheca  matrix.  At  the  start 
of  mating,  the  queen  flies  by  her  maximum  speed.  Then 
a  drone  is  chosen  from  the  population  randomly.  Then 
a  random  number  between  zero  and  one  is  generated. 
If  this  random  number  is  less  than  the  mating  probabil¬ 
ity  evaluated  by  Eq.  (20),  then  the  drone  sperm  is  added 
to  the  queen  spermatheca,  else  another  drone  is  chosen 
randomly  and  the  process  is  repeated. 

Step  8:  Breeding  process  and  updating  drones’  population:  the 
breeding  process  should  be  implemented  as  it  was 
described  in  Section  6.2  by  the  use  of  Eqs.  (22)-(25).  After 
the  termination  of  the  breeding  process,  the  repository 


Fig.  7.  3D  plot  of  the  Pareto-optimal  solutions  found  for  case  study  1  considering  emission,  voltage  deviation  and  cost  objective  functions. 
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Table  7 

Objective  function  values  in  all  cases  (Test  system  1). 


Cases 

Importance 

Wi 

w2 

w3 

w4 

/i 

f2 

h 

U 

Case  I 

- 

- 

- 

- 

85.583 

0.053385 

158.060 

26,215.79 

Case  II 

- 

- 

- 

- 

94.805 

0.048988 

158.429 

25,335.31 

Case  III 

- 

- 

- 

- 

139.533 

0.062180 

154.181 

49,954.62 

Case  IV 

- 

- 

- 

- 

90.596 

0.054528 

158.260 

25,280.76 

CaseV 

0.33 

0.33 

0.33 

0 

101.460 

0.052126 

156.191 

_ 

0.25 

0.25 

0.5 

0 

110.788 

0.054349 

155.402 

- 

0.25 

0.5 

0.25 

0 

92.863 

0.051337 

156.992 

- 

0.5 

0.25 

0.25 

0 

92.863 

0.051337 

156.992 

- 

Case  VI 

0.33 

0 

0.33 

0.33 

85.583 

- 

158.060 

26,215.79 

0.25 

0 

0.25 

0.5 

85.583 

- 

158.060 

26,215.79 

0.25 

0 

0.5 

0.25 

104.620 

- 

156.356 

32,808.73 

0.5 

0 

0.25 

0.25 

85.583 

- 

158.060 

26,215.79 

Case  VII 

0 

0.33 

0.33 

0.33 

_ 

0.050984 

157.698 

27,776.86 

0 

0.25 

0.25 

0.5 

- 

0.050984 

157.698 

27,776.86 

0 

0.25 

0.5 

0.25 

- 

0.055561 

155.870 

35,197.85 

0 

0.5 

0.25 

0.25 

- 

0.050984 

157.698 

27,776.86 

Case  VIII 

0.33 

0.33 

0 

0.33 

94.805 

0.048988 

- 

25,335.31 

0.25 

0.25 

0 

0.5 

94.805 

0.048988 

- 

25,335.31 

0.25 

0.5 

0 

0.25 

94.805 

0.048988 

- 

25,335.31 

0.5 

0.25 

0 

0.25 

85.583 

0.053385 

- 

26,215.79 

Case  IX 

0.25 

0.25 

0.25 

0.25 

94.805 

0.0489880 

158.429 

25,335.31 

0.2 

0.2 

0.2 

0.4 

94.805 

0.0489880 

158.429 

25,335.31 

0.2 

0.2 

0.4 

0.2 

96.610 

0.051103 

156.957 

30,392.10 

0.2 

0.4 

0.2 

0.2 

94.805 

0.048988 

158.429 

25,335.31 

0.4 

0.2 

0.2 

0.2 

85.583 

0.053385 

158.060 

26,215.79 

should  be  updated.  During  each  breeding  for  X*  (the  ith 
drone),  the  position  ofX/  is  updated  as  described  in  Section 
6.2. 

Step  9:  If  all  the  drones  are  checked  go  to  the  next  step,  else  go  to 
step  8. 

Step  10:  Updating  the  repository.  In  this  step  the  repository  is 
checked  so  that  all  the  solutions  stored  in  the  repository 
will  be  non-dominated  solutions. 

Step  1 1 :  Updating  the  queen.  A  new  queen  is  chosen  from  the 
updated  repository  randomly. 

Step  12:  Generating  the  queen  speed.  The  queen  speed  is  generated 
as: 

(28)Squeen  =  rand(.)  x  (Smax  —  ^min)  +  *^min 

Step  13:  Termination  criterion.  If  the  number  of  iterations  reaches 
to  maximum  number  of  the  iteration,  finish  the  algorithm, 
else  go  to  step  6. 


7.  Simulation  results 

In  this  section  the  proposed  method  is  examined  on  two  distri¬ 
bution  networks  as  case  studies. 


Table  8 

Specifications  of  RESs  (Test  system  2). 


Capacity  (kW) 

Type 

Location 

RES  1 

500 

FC 

6 

RES  2 

500 

FC 

53 

RES  3 

500 

Wind  energy 

13 

RES  4 

500 

Wind  energy 

60 

RES  5 

500 

Wind  energy 

18 

RES  6 

500 

PV 

71 

RES  7 

500 

PV 

2 

7 A.  Case  study  1 

The  first  test  system  is  a  12.66kV  system  which  consists  of  32 
buses,  5  sectionalizing  switches  and  5  tie  switches  [15].  The  ini¬ 
tial  power  loss  before  the  reconfiguration  is  202.67  kW.  Table  1 
shows  the  location,  capacity,  and  type  of  RESs.  Also,  the  emission 
factors  related  to  NOx,  C02  and  S02  are  shown  in  Table  2.  The  single 
diagram  of  the  network  is  illustrated  in  Fig.  3. 

Since  MHBMO  algorithm  is  used  in  this  paper  for  the  first  time 
to  solve  DFR  problem,  then  first  of  all,  a  single  objective  opti¬ 
mization  of  the  first  two  objective  functions  (total  active  power 
losses  and  voltage  deviation)  is  evaluated  to  compare  the  perfor¬ 
mance  of  the  MHBMO  algorithm  with  the  other  methods  in  the 
area.  The  results  of  comparison  between  the  performance  of  the 
proposed  method  and  some  famous  evolutionary  algorithms  like 
particle  swarm  optimization  (PSO)  algorithm,  genetic  algorithm 
(GA),  honey  bee  mating  optimization  (HBMO),  etc.,  are  shown  in  the 
following.  The  superiority  of  the  proposed  method  over  the  other 
evolutionary  algorithms  can  be  deduced  easily  from  Tables  3  and  4. 
In  these  tables,  the  single  objective  DFR  is  evaluated  while  the  RESs 
are  neglected.  It  is  obvious  that  the  total  power  losses  and  the 
voltage  deviation  objective  functions  which  are  evaluated  by  the 
proposed  method  are  more  satisfying  than  the  results  of  the  other 
algorithms  and  the  good  ability  of  the  algorithm  can  be  inferred.  In 
Table  5,  the  effect  of  RESs  on  the  DFR  problem  is  investigated.  It  is 
evident  that  the  amount  of  pollution  generated  by  the  total  system 
will  be  decreased  in  the  presence  of  the  RESs.  On  the  other  hand, 
the  use  of  these  sources  of  energy  in  the  system  has  resulted  to  a 
considerable  improvement  of  the  other  three  objective  functions’ 
values.  In  order  to  have  a  more  precise  comparison,  the  average 
value,  standard  deviation  and  the  worst  solution  evaluated  by  the 
different  algorithms  is  shown  in  Table  5.  The  bold  number  in  each 
column  shows  the  optimal  minimum  value  in  that  column.  Com¬ 
paring  the  results  of  Table  5  with  those  of  Tables  3  and  4,  it  could 
be  conducted  that  the  performance  of  the  system  after  using  the 
RESs  has  improved  impressively.  Table  6  shows  the  set  of  non- 


Table  9 

Comparison  of  objective  functions  evaluated  by  different  methods  considering  RESs  for  25  trails  on  Test  system  2. 
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Table  10 

Some  of  the  non-dominated-solution  found  for  MDFR  problem  (Test  system  2). 


Power  losses  (kW) 

Voltage  deviation  (pu) 

Cost  ($) 

Emission  (lb) 

1 

394.845 

0.0518573 

1,155.617 

354,055.00 

2 

411.635 

0.0511236 

1,159.601 

339,811.14 

3 

418.340 

0.0507539 

1,159.301 

342,031.88 

4 

434.698 

0.0518089 

1,156.824 

352,753.40 

5 

470.137 

0.0452303 

1,157.585 

354,907.74 

6 

421.103 

0.0575129 

1,160.041 

338,810.85 

7 

443.167 

0.0574801 

1,161.901 

336,764.20 

8 

423.218 

0.0493687 

1,159.298 

343,783.64 

9 

432.495 

0.0575949 

1,161.843 

336,927.69 

10 

434.013 

0.0518284 

1,157.458 

350,705.02 

11 

484.735 

0.0613938 

1,153.954 

370,825.78 

12 

423.502 

0.0499521 

1,158.755 

345,848.50 

13 

463.612 

0.0597616 

1,155.460 

361,918.89 

14 

428.347 

0.0486935 

1,157.582 

349,087.17 

15 

425.001 

0.0581013 

1,159.491 

341,750.19 

16 

535.687 

0.0597561 

1,167.844 

331,303.69 

17 

421.624 

0.0505263 

1,158.129 

345,749.57 

18 

442.959 

0.0487904 

1,161.286 

338,745.86 

19 

448.565 

0.0501727 

1,155.881 

360,702.48 

20 

428.347 

0.0486935 

1,157.582 

349,087.17 

dominated  solutions  found  for  MDFR  problem  considering  RESs. 
Each  of  these  solutions  can  be  the  best  one,  depending  on  the  pre¬ 
specified  priorities  by  decision  maker.  For  example,  if  the  main  goal 
is  to  minimize  the  power  losses,  so,  the  operator  should  select  the 
best  solution  from  the  first  column  of  Table  6.  It  should  be  noted  in 
this  table,  the  behavior  of  each  objective  function  is  related  to  the 
behaviors  of  the  other  objective  functions.  In  order  to  see  this  rela¬ 
tionship  more  clearly,  along  with  Figs.  4-7,  in  Table  7,  the  results 
of  different  combinations  of  the  objective  functions  are  considered 
as  follows: 

Case  I:  Considering  function /j  ( neglecting  /2,/3  and /4) 

Case  II:  Considering  function  f2  ( neglecting /j,/3  and /4) 

Case  III :  Considering  function/3  ( neglecting /j,/2  and /4) 

Case  IV:  Considering  function /4  ( neglecting /j,/2  and /3) 

Case  V:  Considering  functions /j,/2  and /3  ( neglecting /4) 

Case  VI:  Considering  functions/i,/3  and/4  (neglecting f2) 

Case  VII:  Considering  functions /2,/3  and  f4  (neglecting/!) 

Case  VIII:  Considering  functions /i,/2  and  f4  ( neglecting /3) 

Case  IX:  Considering  functions /i,/2,/3  and  f4 
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As  mentioned  in  Section  5,  a  fuzzy  clustering  technique  is  used  to 
control  the  size  of  the  repository.  Also,  in  this  paper  a  fuzzy  decision 
making  procedure  is  adopted  to  obtain  “most-preferred”  solution. 
In  order  to  aim  this  goal,  the  importance  of  each  objective  function, 
i.e.  w/,  should  be  specified,  such  that  Y^i= iw*  =  1* 

To  have  better  comparison  between  the  single  objective  and 
multiobjective  optimization,  the  results  of  the  single  objective 
optimization  cases  (cases  I— IV)  as  well  as  the  multiobjective  opti¬ 
mization  cases  (cases  V-IX)  are  shown  in  Table  7.  In  the  cases  V-IX, 
the  total  emission,  the  voltage  deviation,  the  total  power  losses  and 
the  total  cost  objective  functions  are  neglected,  respectively.  In  the 
first  part  of  Table  7  (cases  I— IV),  the  value  of  each  objective  func¬ 
tion  after  applying  single  optimization  is  shown.  The  results  show 
that  in  some  cases  the  improvement  of  an  objective  function  will 
affect  the  other  objective  functions  in  the  same  manner  while  it 
may  negatively  influence  another  objective  function.  Some  of  the 
most  important  results  which  can  be  inferred  from  this  table  are  as 
follows: 

-  From  Table  7,  we  can  find  that  in  a  large  range  of  variation, 
the  behavior  of  f2  (voltage  deviation)  and  /4  (emission)  are  the 
same;  that  is  when  f2  is  minimized  individually, /4  is  also  min¬ 
imized  and  when  it  increases  the  other  one  will  increase,  too. 


T.  Niknam  et  al.  /  Journal  of  Power  Sources  196(2011)  8881-8896 


8893 


Switch  3 

Bus  1  /  84 


S/S  I 


A  H-H-t-fc 


JHJ  i 

1 


rc 


I 

I  /  S/S 

I  a  sy  47  i 

H-HH-I-I-H- 


B 


U  |  |  I w  J  |  ,  vvmu 

Vm  .  fiTTr  i 

fS6  ||4  Wind  V  \  *J _ 

To  4)  \»l - 1  JjjN  V  S,  IPV 

I  16  VI  W  /  W  72  bi  . 

I'H1)I|^U  "V-1  kTTi  n  1 1 

U|W1  ^ 

I  -5*  *>n  T4  71 


I  I  «  ^  \  V  I* 

-±j 

'  7+1,  \  117 

HH-1  1 


H-H-H" 

3S  41  42  *3  77 

j  Bus  or  load  center 


43  46 


Sectionalizing  switch 


—  —  Tie  switch 


2 

G 

H 

I 

J 

K 


Fig.  8.  Single  line  diagram  of  85  bus  test  system. 


The  evidence  of  this  claim  is  the  results  of  the  cases  VII,  VIII,  IX. 
In  case  VII,  in  the  2nd  and  4th  row,  when  the  weighting  factor 
of  f2  and/4  is  increased  (so  increasing  their  importance  in  the 
optimization  when  neglecting  /j )  the  value  of  these  two  objec¬ 
tive  functions  is  the  same  (J2  =  0.050984  p.u;  f4  =  27,776.86  lb). 
This  happens  similarly  in  case  VIII,  in  the  2nd  and  3rd  row 


(J2  =  0.048988  p.u;/4  =  25,335.31  lb),  and  in  case  IX,  in  the  2nd  and 
4th  row  C h  =  0.0489880  p.u;  f4  =  25,335.31  lb). 

-  From  this  table,  it  can  be  deduced  that/3  (total  cost)  has  a  con¬ 
flicting  behavior  to  the  other  three  objective  functions  especially 
f2  and  f4.  In  case  VI  (neglecting  voltage  deviation  objective  func¬ 
tion),  increasing  the  importance  of  f3  (by  increasing  w3),  has  a 


Fig.  9.  3D  plot  of  the  Pareto-optimal  solutions  found  for  case  study  2  considering  total  active  power  losses,  voltage  deviation  and  cost  objective  functions. 
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Fig.  10.  3D  plot  of  the  Pareto-optimal  solutions  found  for  case  study  2  considering  total  active  power  losses,  voltage  deviation  and  emission  objective  functions. 


negative  effect  on  the  other  two  objective  functions  (/j  and  /4). 
Similarly  in  case  VII,  the  increasing  importance  of/2  and/4  ( w2  and 
w4)  has  resulted  in  the  same  solution,  while  the  same  increase  in 
the  importance  of  f3  has  resulted  in  a  different  solution  which 
shows  a  conflicting  behavior  with  respect  to  the  other  two  objec¬ 
tive  functions.  The  same  observation  can  be  inferred  from  case 
IX. 

-  The  dependency  of/j  (total  power  losses)  versus  the  other  func¬ 
tions  is  to  some  extent  weak.  From  case  V  it  can  be  figured  out 
that/i  has  a  similar  behavior  with  /2  while  the  confliction  with 
f3  (as  mentioned  before)  is  evident.  In  case  VI,  it  shows  the  same 
behavior  with /4.  But  in  case  VIII,  a  different  conflicting  behavior 
on  both  f2  and  f4  can  be  seen.  Similarly  in  case  IX,  /j  is  in  contrast 


with  the  other  objective  functions.  Therefore,  it  can  be  deduced 
that  /i  behavior  in  progress  is  not  in  direct  relationship  with  f2 
and  /4. 

-  In  case  IX,  3rd  row,  the  main  importance  among  the  four  objective 
functions  is  applied  to  the  cost  function  (/3).  Nevertheless,  the 
value  of  this  objective  function  (cost)  is  more  than  that  of  case  III. 
In  fact  this  difference  in  the  cost  values  of  the  two  cases  indicates 
the  extra  value  which  must  be  paid  as  the  cost  of  decreasing  the 
total  active  power  losses  as  well  as  enhancing  the  security  level 
(voltage  deviation)  and  reducing  the  total  emission. 

-  As  before-mentioned,  each  Pareto  solution  of  the  repository  is 
indicative  of  an  alternative  option  for  the  system  operator  as 
a  decision  maker.  Indeed,  after  applying  all  the  limitations  and 


Fig.  11.  3D  plot  of  the  Pareto-optimal  solutions  found  for  case  study  2  considering  total  active  power  losses,  emission  and  cost  objective  functions. 
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conditions  which  is  determined  by  the  operator  as  the  require¬ 
ments  of  the  system,  it  is  needed  to  extract  the  best  compromised 
solution  according  to  the  preferences  among  the  stored  solu¬ 
tions  of  the  repository.  It  can  be  seen  in  case  IX  of  Table  7  that 
with  the  same  preferences  of  the  objective  functions,  the  best 
compromised  solution  is  94.805  kW,  0.0489880  p.u,  1 58.429$  and 
25,335.31  lb  for/^/2,/3  and/4,  respectively. 

Neglecting/4, /2,/i  and  f3  in  the  multiobjective  optimization,  the 
best  compromised  solution  can  be  achieved  from  the  first  row  of 
the  cases  V  to  VIII,  respectively. 

To  show  the  improvement  of  each  objective  function  with 
respect  to  the  other  objective  functions,  the  three-dimensional 
(3D)  plot  of  the  non-dominated  solutions  of  the  multiobjective  DFR 
problem  are  shown  in  Figs.  4-7.  Each  star  denotes  a  Pareto-optimal 
solution  in  these  figures.  The  dependency  of  each  objective  func¬ 
tion  improvement  to  the  other  objective  functions  can  be  seen  in 
these  figures,  too. 

7.2.  Case  study  2 

The  second  case  study  is  a  standard  11.41<V  radial  distribution 
system  which  is  consisted  of  2  substations,  1 1  feeders  and  85  buses 
and  96  switches  [34].  The  initial  active  power  loss  and  maximum 
voltage  deviation  neglecting  RESs  are  531.99  kW  and  0.052  p.u, 
respectively.  In  Table  8,  the  location,  capacity  and  the  type  of  each 
RES  are  shown.  The  single  diagram  of  the  second  test  system  is 
depicted  in  Fig.  8. 

To  understand  the  behavior  of  each  objective  function  with 
respect  to  the  others,  in  Table  7  the  results  of  a  complete  inves¬ 
tigation  is  shown  and  explained.  So,  in  order  to  avoid  unnecessary 
comparisons  in  case  study  2,  the  performance  of  the  proposed  algo¬ 
rithm  is  assessed  in  the  presence  of  RESs. 

In  Table  9,  the  complete  comparison  between  the  values  of  all 
objective  functions  is  implemented.  All  the  objective  functions  have 


improved  effectively  and  the  satisfying  performance  of  the  pro¬ 
posed  method  is  evident.  Also,  by  comparing  the  values  of  the 
standard  deviation,  mean  value,  time  of  the  run  and  the  worst 
solution  found  by  each  algorithm,  the  superiority  of  the  proposed 
algorithm  is  verified. 

Some  of  the  non-dominated  solutions  found  in  multiobjective 
DFR  problem  considering  the  RESs  are  shown  in  Table  1 0.  In  order  to 
see  the  behavior  of  each  objective  function  according  to  the  others, 
in  Figs.  9-12  the  three-dimensional  (3D)  plot  of  the  non-dominated 
solutions  are  shown.  Each  star  denotes  a  Pareto-optimal  solution. 
The  stars  which  are  shown  by  rectangular  box  in  the  figure  are 
optimal  solutions  which  are  found  during  single  optimization  of 
each  of  the  objective  functions. 

8.  Conclusion 

In  this  paper  an  appropriate  modified  algorithm  based  on  FIBMO 
(MFIBMO)  algorithm  is  introduced  to  solve  the  multiobjective  DFR 
problem  considering  RESs.  The  objective  functions  consist  of  the 
total  active  power  losses,  the  voltage  deviation  of  each  bus  from 
its  nominal  value,  the  total  cost  of  the  system  (including  the  grid 
and  the  RESs)  and  the  total  emission  produced  by  the  system. 
In  the  proposed  method,  a  set  of  non-dominated  solutions  called 
Pareto-optimal  solutions  are  found  and  stored  in  the  repository 
which  its  size  is  controlled  by  the  use  of  fuzzy  clustering  method. 
Also  a  fuzzy-based  decision  maker  was  introduced  to  select  the 
‘best’  compromised  solution  according  to  his/her  experiences  and 
preferences.  To  improve  the  original  FIBMO  algorithm,  the  mating 
process  is  enhanced  so  that  to  improve  the  exploration  ability  of 
the  algorithm  in  the  entire  search  space.  In  order  to  see  the  fea¬ 
sibility  and  ability  of  the  proposed  method,  MHBMO  algorithm  is 
applied  to  two  test  systems  and  the  results  are  compared  by  some 
of  the  most  famous  optimization  algorithms.  The  results  show  the 
good  performance  and  credibility  of  the  proposed  method  in  the 
MDFR  problem.  Also  it  was  shown  that  the  behaviors  of  the  voltage 
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deviation  and  total  emission  objective  functions  are  similar  to  each 
other  while  the  cost  objective  function  has  a  conflicting  behavior 
regarding  the  other  objective  functions.  Also,  the  superiority  of  the 
proposed  method  is  shown  with  respect  to  the  other  optimization 
methods  in  the  area. 
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